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Abstract 

A simple 1-D relativistic model for a diatomic molecule with a double point interaction potential 
is solved exactly in a constant electric field. The Weyl-Titchmarsh-Kodaira method is used to 
evaluate the spectral density function, allowing the correct normalization of continuum states. The 
boundary conditions at the potential wells are evaluated using Colombeau's generalized function 
theory along with charge conjugation invariance and general properties of self-adjoint extensions for 
point-like interactions. The resulting spectral density function exhibits resonances for quasibound 
states which move in the complex energy plane as the model parameters are varied. It is observed 
that for a monotonically increasing interatomic distance, the ground state resonance can either go 
deeper into the negative continuum or can give rise to a sequence of avoided crossings, depending 
on the strength of the potential wells. For sufficiently low electric field strength or small interatomic 
distance, the behavior of resonances is qualitatively similar to non-relativistic results. 
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I. INTRODUCTION 



The interaction of matter with very intense laser field has been a very active field of 
research in the last few decades and led to the discovery of very interesting phenomena such 
as above-threshold ionization, high harmonic generation and others l|, Q]. These studies 
were motivated mostly by the advent of new lasers reaching high intensity levels. In that 
regime, the traditional theoretical tools such as perturbation theory are not applicable and 
thus, many other avenues were explored to study these systems mathematically 2|. More 
recently, the development in laser technologies have allowed to reach unprecedented inten- 
sity levels (10 20 W/cm 2 and higher [3}]) providing an opportunity to study a new range of 
nonperturbative parameters 0]. The interaction of these lasers with matter, such as elec- 
trons, atoms or molecules, allows to probe relativistic effects: the ponderomotive energy of 



a. 



an electron in that regime has an order of magnitude close to the rest mass energy 
In this setting, the Dirac equation should be used to give a consistent description of these 
phenomena instead of the non-relativistic Schrodinger equation. 

From a mathematical point of view, finding a solution to the Dirac equation is a very 
challenging task because of its intricate matrix structure. For this reason, exact solutions 
can be found only for a few special cases describing highly symmetric systems 6^8|. Ap- 
proximate solutions however can be determined through semi-classical techniques and in 
some cases, realistic systems can be analyzed by using these analytical methods Q. An- 
other approach is to use numerical methods; this is the subject of numerous studies for 



both the time-dependent [4| and time-independent cases [10|, lU|. However, these solutions 
often require an important amount of computational resources and the numerical methods 
require special care to circumvent numerical artifacts. Moreover, these calculations are of- 
ten very complicated, albeit being more realistic, and the correct physical interpretation 
is often harder to reach. For these reasons, it may be helpful to consider simpler models 
to understand the basic physical features of a given system. In this article, we adhere to 
this philosophy: we are considering a simplified relativistic model for a diatomic molecule 
in a quasi-static electric field. More precisely, the molecule is modeled by two attracting 
point-like potentials (double delta function potential) and is subjected to a constant electric 
field. This approximates in a crude way the short-range Coulomb potential of the nuclei and 
its interaction with a slowly varying electromagnetic laser field. To simplify the analysis fur- 
ther, we consider a one-dimensional model for which analytical solutions can be determined. 
This last approximation can also be justified in the strong field limit where ionization along 
the electric field is the more important mechanism in which case a "real" linear polyatomic 
molecule behaves almost like a 1-D system [l2[ 13 ]. 

This kind of simplified models using point-like interaction has been studied extensively 
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both relativistically 14j-|l9| and non-relativistically 20|, |2l| for scattering and bound states 



problems when the electric field is absent. In these cases, an analytical solution can be 
calculated easily. In the presence of a constant homogeneous electric field, the mathemat- 
ical problem is more challenging because the differential operator becomes singular on the 
domain limits (at x = ±00), which complicates the boundary condition prescription. Also, 
the continuum spectrum a c of the differential operator changes dramatically: there is a con- 
tinuum on the whole range of energy (a c = K) while bound states become resonances for 
sufficiently small field strength (a resonance corresponds to a pole of the resolvent operator 
on the "unphysical" Riemann sheet 22j). This is to be contrasted with the zero field case 
in which the continuum spectrum is a c = (—00, — mc 2 } U [mc 2 , 00), the bound states in this 
case are points in (—mc 2 , mc 2 ). To cope with these difficulties, analytical methods have been 
developed and applied to these simple models. For instance, Green's function methods were 
used in 
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25] to compute the Stark energy shift and decay rates of quasibound states. The 
Weyl-Titchmarsh-Kodaira (WTK) method 26j, |27|, which allows to compute the spectral 
density of second-order singular operators, was used to analyze many one-dimensional mod- 
els 
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29] and to develop numerical methods for arbitrary 1-D potentials |30|. The same 



kind of approach was used in 3l| to evaluate the non-perturbative Stark effect in hydrogen- 
like atoms. Finally, the 1-D non-relativistic diatomic molecule with point-like interactions 
in a constant field has been investigated in [32j using this WTK technique. In their analysis, 
the authors were able to compute analytically the spectral density and showed results in very 
good qualitative agreement with real diatomic molecules. In our study, a similar approach 
is used to study the latter in the relativistic regime using the extension of the WTK method 
to the Dirac equation [33]. This methodology has been used in [34j-l37j to study various 
potentials, while the atomic case with a single delta function potential was treated in 38]. 

An important advantage of working with point-like potentials is that they can be char- 
acterized by boundary conditions on the Dirac potential well positions. This procedure in 
the non-relativistic case is well-known and a rigorous discussion of this subject can be found 
in [21]. For the Dirac equation, this has been the subject of debates in the last few decades 



because of ambiguous results: different methods yielded different boundary conditions [16 ]. 
This ambiguity was related to the fact that the wave function has a jump discontinuity and 



is multiplied by a delta function, resulting in an ill-defined product of distributions [39|, |40|. 
This was also studied rigorously by looking at the self-adjoint extensions of the Dirac op- 
erator and it was shown that there exists a four-parameter family of self-adjoint operators 
corresponding to point interactions [l?], 0-44]. In this work, we make use of certain results 
rom the algebra of Colombeau's generalized function to deal with products of distributions 
45 48|. By combining these ideas with general properties of self-adjoint extensions 41] 



and charge conjugation invariance, we show that it is possible to select a specific boundary 
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condition which corresponds to an electrostatic Dirac potential wells. 

This article is separated as follows. In Section [Til the WTK method is shortly presented. 
This discussion focuses on general results which allow to comprehend the general ideas and 
to understand how to apply the methods in our case. In Section IHIl the free case (with 
no electric field nor point-like potential) is considered to show an implementation of the 
WTK method and an important qualitative difference with the non-relativistic case, i.e. 
the appearance of the negative energy continuum. The diatomic molecule without electric 
field is considered in Section [TV] as a consistency check for the full case, which is treated in 
Section IIVI The calculation of boundary conditions using Colombeau's theory is relegated 
to Appendix along with a few comments on the Dal Maso-Le Floch-Murat definition of 
products of distribution as a bounded Borel measure. Note that throughout this work, we 
are working in units where c = m = h = 1, so that the Compton radius is unity (m and c 
are kept explicitly in all equations to easily switch to atomic or natural units). 



II. GENERAL THEORY OF EIGENFUNCTION EXPANSION FOR THE DIRAC 
EQUATION 



In this section, a brief review of the WTK method is presented. This discussion focuses 
on main results and we refer the interested reader to the extensive mathematical literature 
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49 



50] 



on the subject for more details 

To begin, let us consider an equation defined on a domain x G M of the form 

Htfj(x) = Eip(x) (1) 

where H is the Dirac operator, E is the eigenenergy and V G L 2 (R)®C 2 is the wave function 
(bi-spinor). The main objective here is to expand the wave function over an eigenfunction 
{4>e)e basis. Assuming that if is a self-adjoint linear operator on T>(H) defined in a Hilbert 
space, the spectral theorem gives the spectral decomposition: 



H= / EdP E (2) 

where u(H) is the spectrum of operator H and P% is a projection operator onto the 
eigenspace spanned by the eigenvectors ipE- Written in this abstract form, the spectral 
theorem is of limited utility for practical applications: its mathematical properties can ob- 
viously be studied but it says little concerning the explicit expression of Pe- 

The WTK method provides an explicit realization of the spectral theorem for the Dirac 
operator in 1-D: it allows, from the knowledge of a solution of the equation, to construct a 



generalized eigenfunction expansion, in the form of ([2]) 5l|. At the same time, it yields the 
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spectral density which has a very important physical meaning since the latter exhibits the 
main features of the spectrum (bound states, resonances and continuum part). Finally, the 
WTK method can be used to treat "singular" problems where the operator H has a singu- 
larity on one or both domain boundaries. This case is relevant when the system is immersed 
in a static electric field, as considered in this study, and where we have lim^i V(x) = =Foo 
for the potential V. The main result of the WTK analysis is that any smooth function / 
can be expanded as 

/oo 1 
J2li(^E)p ij {E)f j (E)dE (3) 

where Pij(E) are the components of the matrix-valued spectral density, j (x,E) := ■u ± (x), 
7i(x, E) := v ± (x) (the functions u and v will be defined later, see (|T2|) and (jTB"]) ) and 

/oo 
jj(x,E)f(x)dx (4) 
-00 

are Fourier-like coefficients. Note here that by applying the Hamiltonian H to /, one obtains 
an expression for the projection operator P% in terms of the spectral density. The latter can 
be specified explicitly because they can be related to the classical solution. This procedure 
will be discussed further in the following. 

The WTK method was originally developed to extend the Sturm-Liouville problem to 
infinite domain with singular boundary conditions, but it can be carried to the relativistic 
case and more generally, for the solution of any first order systems of equation of the form 
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p{x)ip 1 {x) + q(x)ip 2 {x) + d x ip2(x) = Eipi(x), (5) 
q(x)ipi(x) + r(x)4j 2 {x) - d x tpi(x) = Eip 2 {x), (6) 

for all x G M, where p, r, q are real functions which can have singularities on the domain 
boundaries. The 1-D Dirac equation is a special case of the last equation. More specifically, 
we first consider the case where x~ < x < x + (we will let x~ — > —00 and x + — > 00 later) . 
The Dirac equation is given by 

Eip(x) = [—icad x + f3mc 2 + V(x)] i[>(x), x G (7) 

where in 1-D, ip is a bi-spinor, V is the scalar potential, m is the fermion mass and E is the 
energy. It is convenient to work in a representation where the Dirac matrices are given by 

a = —o y and (3 = —a z (8) 
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which obey the appropriate anticommutation relation^]. This representation ensures that 
we start with an equation of the form of fl5J) and with real coefficients. Explicitly, we 
have for all x G [x~, x + ] 

E 1 

— ipi(x) = d x ip2(x) — mcipi(x) H — V(x)ipi(x), (10) 
E 1 

—ijj 2 (x) = —dxtpxix) + mc^ 2 (x) + -V(x)^ 2 (x). (11) 
c c 

Now, let u + ,v + and u~,v~ defined for x G [x ,x + ] and x G [x~,x ] respectively (for x G 
[x~,x + ]), be two solutions of ( Til)]) and (flTl) having the following "initial data" [331] 

<(x ) = l, 4(^o) = 0, (12) 

uf(x o ) = 0, ^(x ) = l. (13) 

Here, xq can be chosen arbitrarily and from now on, we set xq = 0. This particular set of 
"initial" data guarantees that the condition u^(x)vf(x) — uf(x)vf(x) = 1 is fulfilled for all 
values of x, respectively (it can be shown that the value of this expression does not depend 
on the space coordinate) and that the solutions u and v are linearly independent. Thus, 
a set of two eigenf unctions is found with initial conditions (at an arbitrary point within the 
domain) chosen such that the Wronskian is normalized to 1. Thus, we can write the general 
solution as a linear combination of these solutions (up to an overall normalization constant) 
such as 

ijj+(x) = u + (x) + l + (E)v + (x), xe[0,x + ], (14) 

i/j~(x) = u~{x) + r(E)v~(x), xe[x~,0], (15) 

where l^ 1 G C are integration constants that need to be fixed by boundary conditions at 
x = x , respectively. The choice of these boundary conditions is very important to ensure 
the self-adjointness of the operator H. A general choice is given by 

V'i : (z ± )coB/3 + Vj(ar ± )siii/3 = (16) 

where (3 G M is an arbitrary parameter. This choice is very similar to the one in Sturm- 
Liouville problems: it contains the Dirichlet and Neumann conditions as special cases, and 
it ensures that the operator is self-adjoint. The integration constants obey ^(E) x ~^ ±oc ) 
m ± (E). The value of m ± (£') such that (fj"6"|) are fulfilled lies on a circle in the complex plane 
parametrized by (3 and two outcomes are possible 26|, |49) : 



1 This can be easily obtained from the usual Dirac representation by using a change of representation 
transformation defined as 

U= -j=(a x - a y ). (9) 

such that a ncw = f/aairac^ j Aicw = C^/?dirac^^ and UU' f = I. Recall that the Dirac equation is invariant 
under these transformations. 
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1. The radius of the circle vanishes (limit-point). 



2. The radius of the circle stays finite (limit-circle). 



In the former, it can be proved (WTK theorem) that there exists a unique non-trivial 
solution (independent of (3) such that ^ ± G L 2 (0, oo) <g> C 2 (L 2 (— oo, 0) <8> C 2 , respectively) if 
Q(E) > 26|, |49|. In the following, we consider only this limit-point case since the system 
under study falls into this category [33|. This can be seen easily by looking at the explicit 
solutions ^ calculated in the next sections and by noting that ip + (jL L 2 (0, oo) (g) C 2 and 
%p~ (jL L 2 (— oo, 0) <8> C 2 when $s(E) = 0. This, by definition, is the limit-point type 49]. 

The functions m ± (E) are the Weyl-Titchmarsh m-functions and their knowledge allow us 
to compute the spectral density. The general procedure to compute m ± (E) is to construct 
a linear combination of two solutions obeying (I12p and (1131) . Then, the functions m ± (E) 
are chosen such that the linear combination vanishes at the boundaries when the energy has 
a non-zero positive imaginary part. For example, if the singularity is at x = oo, one must 
look at the asymptotic solution and make sure it vanishes in that limit when ^s(E) > 0, 
guaranteeing that the solution is in L 2 (0, oo) ® C 2 and obeying the WTK theorem. 

The spectral density described earlier in ([3]) can be related to the Weyl-Titchmarsh 
functions by applying the the resolvent operator on / and by comparing with usual results 
for the Green's function {|3oj]. The final result is that the spectral density is given by 



P{E) 

where the matrix M is defined as 



M(E) 



1 



m+(E) -mr(E) 



lim -QM(E + ie) 

e^0+ 71 



-1 



| [rn+(E) +mr(E)] 
\ [m + (E) + m~(E)] m + (E)m-(E) 



(17) 



(18) 



However, it is more convenient for our purpose to work with the "trace spectral density" 



given by p = Tr(p), or more explicitly by 30_|, |32] 



p(E) := lim — S 

e^0+ 71 



m + (E + ie)m~(E + ie) + 1 
m + (E + ie) — m~(E + ie) 



(19) 



The latter contains all the physical information on the system, i.e. the whole spectrum is 
included in this expression 49_|. Note here that if the potential is an even function, the 



relation m + (E) = —m (E) holds and we have 

1 



p(E) 



lim — S 

>o+ 2tt 



m + (E + ie) + 



1 



m + (E + ie 

Thus, for this specific case, we only have to construct the solution ip 



(20) 
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III. FREE CASE 



Now, we consider an explicit example of the theory described in the preceding section: 



the free case on the rea. 
can be found in [l?], |41| 



line. This is probably the simplest possible case and the solution 
although in the form of the resolvent operator. The reason for 
presenting this example is twofold: it allows to implement the method in a simple setting 
and also, it will serve as a validation tool for the results obtained in the next sections. 
Specifically, the free Dirac equation is given by 

E 

— ij) x (x) = d x i) 2 {x) - mc^i(x), (21) 

—ip 2 (x) = -d x ipi(x) + map 2 (x). (22) 

The general solution to this system of equation can be easily computed in terms of trigono- 
metric functions. This solution, for all x, is given by 

ipi(x) = c\ sin (^-xj + c 2 cos (j-x^j , (23) 

ifa(x) = -ci— — — =■ cos (-x) + c 2Tr -^ — ? sin ( -x ) , (24) 
E — mc z V c / E — mc z V c / 

where p = \JE 2 — m 2 c 4 (note that we have p = iy/m 2 c i — E 2 when \E\ < mc 2 , that is 
we chose the Riemann sheet where ^s{p) > 0) and are integration constants. These 
integration constants have to be fixed by using suitable boundary conditions. According to 
the theory developed in the last section, we are interested in finding solutions obeying the 
boundary conditions in (1T2"]) and (IT51) . We define these solutions as u and v, and they are 
given explicitly for all x by 

Ui(x) = cos (J- X ^j ) (25) 

u 2 {x) = —— — - sin (~x) , (26) 
E — mc z V c / 

V\[x) = sin y-xj , (27) 

^2^) = cos (-XJ . (28) 



.c 

From these solutions, we now form the general solution ip + (in the free case, the potential 
is even, so ip~ is not needed). It is given by 

i/j+(x) = u(x) +m + (E)v(x), xeR + . (29) 

For general m + however, this solution is not in L 2 (0, oo) ® C 2 when ^s(E) > 0, as can be seen 
by looking at the function behavior at infinity: asymptotically, we have tp + (x) ~ e ±l ^ x . 



s 



However, by following the WTK prescription and giving a small positive imaginary part e 
to the eigenenergy, we have 

p = y/(E + ie) 2 - m 2 c 4 = VE 2 - m 2 c 4 + % - = =e + Q(e 2 ) (30) 



yE 2 — iii-c ' 



and thus, there are three possible cases: 
1. If £7 > mc 2 , then 



x— >+oo x— >+oo 



lim e*c x = and lim e ^ x = oo (31) 

2. If E < -mc 2 , then 



lim e i " x = oo and lim e~^ x = (32) 

x— >+oo x— >+oo 



3. If —mc 2 < E < mc 2 , then 



lim e^ x = and lim e"^ x = oo (33) 

x— >+oo x— >+oo 

Ensuring that the diverging exponentials are canceled can be done by choosing the Weyl- 
Titchmarsh m-function to be 

a-/ s I —i tp P 2 , for E > mc 2 , , 



v 



for < —mc 2 



' E—mc 2 ' 

The spectral density of the free case can then be easily computed from (I2"0"j) and we get 

, . fi^, for I El >mc 2 

P(E)= Y f J" 2 " ( 35 ) 

I 0, for \E\ < mc 

This spectral density is plotted in Fig. [TJ Note the appearance of the negative energy 
continuum for E < —mc 2 = — 1. 



IV. DOUBLE DELTA POTENTIAL WITHOUT ELECTRIC FIELD 

The next case considered is the one of a double-delta potential without electric field. This 
represents physically a simple model for the ground state (symmetric) and first excited state 
(anti-symmetric) of a diatomic molecule. We recall to the reader that the same system was 



treated non-relativistically in 32]. In this case, the spectrum has a continuum component 
for E > while bound states appear as Dirac delta peaks at the bound state energies. 
Similar features are found in the relativistic case. The main differences however are the 
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-0.5 



FIG. 1: Free spectral density. 



appearance of the negative energy continuum while the bound states are positioned in the 
mass gap (—mc 2 ,mc 2 ). 

Specifically, the main goal is to compute the spectral density for a potential given by 



V(x) = -g [S(x -R) + S(x + R)] , x e 



(36) 



where 2R is the interatomic distance and g is the strength of the Dirac delta potential wells. 

The calculation proceeds as in the free case: the first step is to find two solutions obeying 
the boundary conditions f|T2|) and (JIB"]) , then one constructs a general solution in L 2 (0, oo) <S> 
C 2 from which we can obtain the m-functions and thus, the spectral density. 

The main difference with the free case is the treatment of the point interaction. As usual 
with these kinds of potential, the solution can be found by solving the free Dirac equation 
by including boundary conditions at x — ±R. However, it is well-known that there is an 
ambiguity in the definition of the boundary condition corresponding to the potential in (|36| 
in the Dirac equation 16J|. This ambiguity can be traced to the fact that the wave function 
should have a jump continuity at x = ±R and thus, the terms ip(x)5(x±R) appearing in the 
equation are actually products of distributions (or generalized function). The appearance of 
the jump discontinuit y in the wave function can be understood heuristically in the following 
way, as suggested in [39|. Let us consider a massless Dirac equation at E — (a finite E 
does not change the argument) and only one Dirac delta potential positioned at x = such 
that —iad x if}(x) = 5(x)ip(x). For this equality to hold, the wave function has to be of the 
form ij; ~ H (where H is the Heaviside function) which implies a jump discontinuity. 
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The strategy used in this work to obtain the boundary conditions is to formulate every- 
thing in the weak sense and consider the product of the potential and wave function as a 
product of distributions from the outset. The main problem is that the set of distributions 
is not an algebra 52j. There exists however a mathematical construction that includes these 



45M48]. 



product of distributions rigorously: the Colombeau's theory of generalized function 
In that framework, the value of a distribution product depends on the "local" structure 
of the distributions considered, i.e. the choice of the limiting function. Using these ideas, 
we derive in Appendix |A] a family of boundary conditions consistent with (IA18j) and thus, 
having the right mathematical properties (self-adjointness and non-relativistic limit). From 
this analysis, the following boundary conditions are obtained: 

,2 



and 





9 2 
1 + — 
4c 2 _ 


-1 r 


M±R + ) = 


g 2 ~ 
1 + — 

. 4c2 _ 


-1 r 


M±R~) = 


r a 2 ~ 
1 + — 

4c 2 _ 


-1 - 


M±R~) = 


\ a 2 ' 
1 + — 

4c 2 


-1 r 



1 - 



9 
4c 2 

l_ 

4c 2 



9 

Ac 2 

l_ 

Ac 2 



M±R~ 

M±R" 
M±R" 



- 9 -^{±R~ 

c 
c 



+ -M±R" 

c 

- ^M±R' 

c 



(37) 
(38) 

(39) 
(40) 



where we defined ip(±R + ) — h m e-s>o if>{±R + e) and ijj(±R~) = lim e ^ ?/;(±_R — e). 

Having now all the necessary elements, it is possible to compute the spectral density for 
the double Dirac system without electric field. The procedure is now summarized: 

1. Compute u and v for x G [0, R), which are solutions obeying the boundary conditions 
(IAT8D and (jA"l9j) ). 

2. Compute m > and v > for x G (R,oo) using the conditions (1371) and (|38l) . This yields 
the solutions on x G [0, 00) as (see Figure [2] for the domain of each function) 



u + (x) 
v + (x) 



u(x) for < x < R 

m > (x) for R < x < 00 

v(x) for < x < R 

v > (x) for R < x < 00 



(41) 
(42) 



3. Compute m + by ensuring that ip + (x) = u + (x) + m + v + (x) is in L 2 (0, 00) <g> C 2 when 
$f(E) > 0. 



4. Compute p(E) from m + . 
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u{x) 
v"(x) 

-R 


u\x) 
v\x) 

R 


it(x) 
v\x) 




w(x) 
v(x) 


m1(x) 
v\x) 



FIG. 2: Domain of the solutions. 



A. Solution in terms of trigonometric functions 



The first part of the calculation was already performed in the free case. In the latter, 
the free Dirac equation, ( 121 j) and (I22p . was solved with suitable boundary conditions. The 
solution obtained in (1251) to (|28|) thus corresponds to the needed solution in the region 
x G [0,R), i.e. w(x) and t>(x). 

The next step is to compute the wave functions in the region x > R. They are given by 



u%(x) 

v i( x ) 
v%(x) 



ai sin I -x I + a 2 cos I —x ; . 
-ax— -cos -x + a 2 — ^sm -x , , 



6i sin ^-^j + b 2 cos ^-x 



.c 

&i— — — - cos ( -x) + 6 2 

— mc z \c / 



sin -x 

£/ — mr V c 



(43) 
(44) 
(45) 
(46) 



The value of the integration constants a 1; a 2 , b± and b 2 can be determined from the boundary 
conditions in ( 137|) and ( |38|) . These conditions yield the following relations: 



u>(R) 
u>{R) 
v>(R) 
v>(R) 



1 + i- 

4c 2 



-i 



#1 

4c 2 

4c 2 

g 2 ' 
i + — 

4c 2 



1 + 



1 + 



9 

Ac 2 

l_ 

Ac 2 

l_ 

Ac 2 

l_ 

Ac 2 



ui(R) - -u 2 (R) 
c 

u 2 (R) + -u^R) 



vi(R) - -v 2 (R) 

c 



v 2 (R) + ^ Vl (R) 



:= G: 



G, 



G v ■ 



(47) 
(48) 
(49) 
(50) 
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As usual, it is possible to solve for the coefficients and we get 




(51) 



(54) 



(53) 



(52) 



This completes the derivation of the solution on x G [0, oo) for u + and v + . The solution on 
x < could be calculated in a similar way but can be easily obtained from ip + owing to the 
symmetry of the potential. The next step is the evaluation of the m-function. 

B. Evaluation of the Weyl-Titchmarsh m-function 

In this section, the function m + is evaluated. Remember that the general solution ip + = 
u + + m + (E)v + should be in L 2 (0, oo) Cg>C 2 when $s(E) > 0. A similar calculation was carried 
in the free case so the details are not repeated here. The final result is that 



This choice ensures that the diverging exponentials are canceled as x — > oo and $s(E) > 0. 
Thus, by combining this last equation with the solution computed in the last section, we 
can evaluate the spectral density. 

C. Evaluation of the spectral density 

The spectral density is given by (1201) when the potential is even. Having determined 
the m-function in the last section, it is now possible to evaluate the spectral density. It is 
a complicated function of the interatomic distance and potential strength, so a numerical 
evaluation seems better suited than showing the equation. However, the bound state energies 
have relatively simple expressions: they are given by the position of the poles in the spectral 
density and thus, satisfy the following transcendental equation 



for the energy E having a value in the mass gap, that is E G [—mc 2 ,mc 2 }. Here, we 
defined p = y/m 2 c A — E 2 and the cases ± correspond to the ground and first excited states 




(55) 




(56) 
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respectively. It was verified that the solutions of (|56|) corresponds to the position of the 
bound state in Figs. EJandlH 

In Fig. [3j the spectral density is plotted for different values of the potential strength 
(g) and a fixed atomic interatomic distance (R = 1). As shown in the figure, there is only 
one bound state at low potential strength (the ground state). When the potential strength 
reaches g = —Amc 3 R + 2cy/Am 2 c 4 R 2 + 1, the second bound state appears (the first excited 
state), while the ground state has a lower energy. By increasing g furthermore, the two 
bound states have lower eigenenergies until the ground state disappears in the negative 
energy continuum when g = Amc'R + 2c\/4m 2 c 4 R 2 + 1. To summarize, we have 

1. For < g < —Amc'R + 2cV4m 2 c 4 i? 2 + 1, the ground state is the only bound state in 
the mass gap. 



2. For -Amc'R + 2cVAm 2 c 4 R 2 + 1 < g < Amc'R + 2cVAm 2 c 4 R 2 + 1, both the ground 
and the excited states are bound states in the mass gap. 



3. For Amc 3 R + 2c\/Am 2 c 4: R 2 + 1 < g, the excited state is the only bound state in the 
mass gap, the ground state has vanished in the negative energy continuum. 

Note here that for g ~ or g ^> 1, the free spectral density is recovered. 

In Fig. HI the spectral density is shown for many values of interatomic distance (R) and 
for fixed value of potential strength (g = 2.0). For small value of R, the energy difference 
between the ground and excited states is larger. As the interatomic distance increases, the 
ground state eigenenergy also increases while the excited state eigenenergy diminishes until 
they merge and form a quasi-degenerate state (at large R). Note also the qualitative change 
in the positive and negative continua: as R gets larger, the oscillating behavior is amplified. 
This phenomenon was interpreted physically in 32J as the resonant backscattering between 
the two potential wells (the so-called Ramsauer-Townsend resonances). 



V. DOUBLE DELTA POTENTIAL WITH ELECTRIC FIELD 



In this section, we consider a more interesting system: the double delta function with an 
electric field. This system can be seen as a very simplified model of a diatomic molecule in 
a slowly varying electromagnetic field. Specifically, the main goal is to compute the spectral 
density for a potential given by 

V(x) = -g [5(x -R) + 5(x + R)] -Fx, iGK (57) 

where 2R is the interatomic distance, F is the electric field and g is the strength of the 
Dirac delta potential wells. Clearly, this potential is singular at x = ±oo and therefore, the 
general theory of eigenfunction expansion should be used. 
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(e) (f) 

FIG. 3: Spectral density for a two Dirac delta potential for varying values of g. The value 
of the interatomic distance is fixed to R — 1.0 while the potential strength is (a) g = 0.2, 
(b) g = 0.5, (c) g = 1.0, (d) g = 3.0, (e) g = 10.0 and (f) g = 50.0. 
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(a) (b) 
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(c) (d) 

FIG. 4: Spectral density for a two Dirac delta potential for varying values of R. The value 
of the potential strength is fixed to g — 2.0 while the interatomic distance is (a) R = 0.1, 

(b) R = 0.5, (c) R = 1.0 and (d) R = 2.5. 

As described above, the first step is to find two solutions obeying the boundary conditions 
( TT2|) and ( TT3|) . To do this, we start by finding the general solution in region —R < x < R and 
x < —R and x > R, i.e. when there is no delta wells. Once we know the general solution, 
the case with delta wells can be easily treated by using suitable boundary conditions at 
x = ±R. The boundary conditions on the point interaction are described in Appendix |A] 



A. Solution in terms of Parabolic Cylinder functions 



The general solution in a constant electric field is known, it was solved in 33], but we 
recapitulate the main steps here for completion. A similar calculation was also performed 
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in 38[. Following 33|, we start by rewriting the spinor components as 
such that 



-j= (ipi + itp 2 ) , ^2 = -j= (ipi - #2) 



(58) 



(59) 



Substituting the latter in f[10|) and ffTTj) . adding and subtracting the two equations, we find 
that 



V>i = 4= (^l + ^2) , ^2 = -4| (fii - fi 2 ) , 



E - V(x) 



E - V(x) 



fli(x) — id x Qi(x) — mcfl 2 (x) = 0, 
Q 2 (x) + id x Q 2 (x) — mcQi(x) = 0. 



Then, solving for Q 2 in (160]) and substituting in (151]) . we get, after a few manipulations 



SIM 



mc 



+ id x 



E + Fx 



F 

■ 11 
1 m c 



fii(x) = 0. 



(60) 
(61) 

s 

(62) 
(63) 



The second equation can be simplified further by letting y = e l * { E+ f X ) > an< ^ this 
yields 



-y + a 



fii(y) = o, 



(64) 



where a 



i^jpr- — \. The last equation has well-known solutions and can involve any 
pairs of independent solutions a mong U(a, ±z),V(a, ±z) or U(—a, ±iz) where U and V are 
the parabolic cylinder functions 53 



54| . However, it is convenient to use a "numerically 
satisfactory pair of solution", that is a pair for which one solution is dominant while the 
other is recessive as — > ±00. This facilitates the calculation of the m-functions. In our 
case, the phase of the argument is arg y = —tt/A and thus, a numerically satisfactory pair 
of solution is given by 54 ] 



fii(y) = C\U{a,y) + c 2 U(-a, -iy), 



(65) 

where U is the parabolic cylinder functions defined in [53J] and Ci t2 are integration constants. 
From this, we can evaluate the second component. In the y coordinates, it is given by 



I 2F 

mc V c 



(66) 
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Then, using the recurrence relations for parabolic cylinder functions 53] 

^zU{a,z)+d z U(a,z) = - ( a +^J U{a+l,z), (67) 
^zU(a,z)-d z U(a,z) = U(a-l,z), (68) 

we get that 



Q 2 (y) = mc x l^-e l '^ Cl U{a + l,y) + — J—e-^c 2 U(-a - 1, -iy). (69) 
lb mc V c 



Here and in the following , we define 







:= U(a 


y), 




u 2 


[a,y) 


:=U(- 


a, -iy), 




Ui 




:= mc^ 




+ 1,2/), 


u 2 








-a - 1,2/), 



(70) 
(71) 

(72) 
(73) 



which allows us to write the general solution as 



^i(x) = ciUi[a,y(x)) + ciUx[a,y(x)) + c 2 U 2 [a,y(x)) + c 2 U 2 [a,y(x)), (74) 
itp 2 (x) = ciU![a,y(x)] - ctU^a, y(x)) + c 2 U 2 [a,y(x)] - c 2 U 2 [a,y(x)], (75) 

(the factors l/v2 were absorbed by redefining integration constants). Following the same 
procedure as in the preceding cases, the integration constants have to be determined from 
boundary conditions (( fT2|) and ( fTBl ) yielding u and v for x G (—R,R). The conditions in 
f )37|) and fl38|) will serve to find m > and i> > for x G (-R, oo), while f )39|) and (|4T)l) will serve to 
find m < and f < for x G (— oo, — i?). 

It is now possible to evaluate the integration constants such that the boundary conditions 
( |T2l) and ( |T3|) are obeyed. Letting j/ : = e - *^ ) be the value where x = 0, we find the 

coefficients 

1 U 2 (a,y ) -U 2 (a,y ) , . 
C\ = —-— ~ , (76) 

2 U 2 {a,y )U 1 {a,y ) - U 2 (a,y )Ui(a,y ) 

C2 = 1 Ui(a,yo) ~ Ui(a,y ) ^ 

1U 2 (a,y G )Ui{a,y Q ) - U 2 (a, y )Ui{a, y ) ' 



and 



c / = £ U2(a,yo) + U^jhyo) ^ 

' 2 U 2 {a, y Q )U x {a, y ) - U 2 (a, y )Ui(a, y ) ' 

c / = _£ Ui(a,yo) + Ui(a,yo) ^ 

2f/ 2 ( a ,?/ )L r i(a,?/o) - ^2(a,2/o)f^i(a,2/o)' 
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for the first (u has integration constants Ci^) and second (v has integration constants 
c[ 2 ) boundary conditions respectively. The next step is to compute the solution for x G 
(—00, —R) and x G (R, 00). 

The general solutions for x G (R, 00) can be written as 



u i( x ) 

Vi(x) 



iv%(x) 



bxU 1 [a,y{x)\ + bxUx[a,y(x)] + b 2 U 2 [a,y(x)) + b 2 U 2 [a,y(x)], 
bxU 1 [a,y{x)\ - bxU 1 [a : y{x)\ + b 2 U 2 [a,y(x)] - b 2 U 2 [a,y(x)], 
biUx^yix^+b^Uxi^yix^ + b^U^yix^ + b^U^yix)], 
biUi[a,y(x)} - tfjj x [a,y(x)\ + b' 2 U 2 [a,y(x)) - b' 2 U 2 [a,y(x)), 



(80) 
(81) 
(82) 
(83) 



while for x G (—00, —R), we have 



uf(x) 
iu 2 (x) 

vf(x) 
ivo(x) 



a 1 U x [a,y{x)\ + axUx[a,y{x)\ + a 2 U 2 [a, y{x)\ + a 2 U 2 [a, y(x)], 
a 1 U x [a,y{x)\ - a\U\[a,y{x)\ + a 2 U 2 [a,y(x)] - a 2 U 2 [a,y(x)}, 
a' x U\[a,y(x)\ + a! x U\[a,y(x)\ + a' 2 U 2 [a,y(x)\ + a' 2 U 2 [a } y(x)}, 
oi?7i[o,y(a:)] - a' 1 t/i[o,i/(a:)] + yfa)] - a' 2 U 2 [a,y(x)], 



(84) 
(85) 
(86) 
(87) 



where the new coefficients a 12 , a[ 2 , 6^2, 6^ 2 can be determined from the boundary conditions 
in (I37p to (140]) . First, let us define the following constants: 



G v :- 



Hi-.- 



9 

Ac 2 

l_ 

Ac 2 
,2 



-1 



1 + 



1 + 



9 
Ac 2 

l_ 

Ac 2 



-1 



-1 

l-^u 2 ^R)± 9 -u i2 (B) 

1 - B) v ^ R) ± 9 c v ^ {R) 

l ~A^) U ^- R ) ±9 - U ^- R ) 



l-^)v 1 , 2 (-R)±^v 2)1 (-R) 



(88) 
(89) 
(90) 
(91) 



The boundary conditions at x — R now become 



u>(R) = G- , u>(R) = G+, 
v>(R) = G~ , w>(fl) = G+ 



(92) 
(93) 



while at x — — R, they are given by 



= . u<(-R) = H u 



v<(-R) = Hl 



-R) = H- 



(94) 
(95) 
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It is now possible to solve for 01,2, a[ 2 , b± j2 , b[ 2 from these equations. We get 

b _ 1 iG+U+[a,y(R)]- g^TMg] 

2 2 ^[a, y(fl)]tf 2 [a, - E^a, y(R)]U 2 [a, y(R)} ' 

fe , _ 1 iG+U+[a,y(R)] ~ G-U^[a,y(R)} 

2 2 U^a, y(R)]U 2 [a, y{R)} - U^a, y{R))U 2 [a, y{R)} ' 



(96) 
(97) 



a 1 iH^maM-Kjl -H+U 2 -[a,y(-R)} 

"~ 2^ 1 [ a , 2/ (- J R)]E/ 2 [a,3/(- J R)]-E/i[a,i/(- J R)]^ 2 [a,y(- J R)]' 

fl / = 1 iH-U+[a,y(-R)]-H+U2[a,y(-R)] m 

1 2^ 1 [ a ,y(-i2)]^ 2 [fl^(-i2)]-^ 1 [a^(-i2)]C)- 2 [a,y(-i2)]' 

where U^ 2 '■= ^1,2 ± ^1,2 and similar expressions for bi,b[, a 2 , a' 2 . However, as will be seen 
in the next section, the latter are not useful for the spectral density calculation. 

The preceding calculation allowed us to evaluate the solutions and v^. Remember 
that these solutions are defined on x > and x < (see Figure [2] for the domain of each 
function). Moreover, these solutions still obey the boundary conditions of the WTK method, 
that is ffT2l) and ffT3l). and can be written as 



u{x) forxG[0,.R) J u{x) for x G 0] 

IT (x) = <( _ , U (X) — < _ , (100) 

u > (x) for x G (R,oa) Im < (x) for x G (— 00, — R) 

Ju(x) forxG[0,i?) Jv(ar) for x G (-#, 0] 

■ir (x) = < _ , v (x) = < . (101) 

lw > (x) for x G (-R, 00) lt> < (x) for x G (— 00, — R) 

We now have all the ingredients to compute the Titchmarsh m-functions and the spectral 
density. 



B. Evaluation of the Weyl-Titchmarsh m-functions 

The m-functions are constructed such that the solutions ip = + m ± (E)v ± are finite 
when x — > ±00 (with the energy obeying ^s(E) > 0). This can be evaluated by looking 
at the asymptotic expansion of the solutions. Asymptotically, in the limit z ^> \a\, the 



parabolic cylinder functions have the following behavior 53J: 



11 3 
U(a,z) ~ e"5 2 V a -^ for |arg(z)| < -vr, (102) 

U(a,z) ~ e -\ z \~ a -\ ±i^^—e^ a e\ z2 z a -\ 

r( 5 -o) 
1 5 

for -7T < ±arg(z) < -7r, (103) 

where T is the usual Gamma function. From these relations and the fact that y = 
e-*5«/f (^p), we find that for the solution ^+ to be in L 2 (0, 00) ® C 2 when 3(E) > 0, 
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it should be 

lim if>+(x) = lim (61 + m + b[) \ Ufayix)] + U l [a,y{x)} \ = 0, (104) 



lim if>£(x) = -i lim (61 +rn + b[) \Ui[a,y(x)] + U 1 [a,y(x))\ = 0, (105) 

as x — > 00, whereas for the solution ip~ to be in L 2 (— 00, 0) ® C 2 when > 0, it should 

be 



lim %p 1 (x)= lim (a 2 + m a' 2 ) < C^fa, y(a;)] + t^K y{x)] \ = 0, (106) 
lim ip 2 ( x ) = ~ i nm ( a 2 + fn~ (^2) \ ^[o, + ^[a, 2/(2)] f — 0, (107) 



as x -)■ — 00. The other terms in the solutions would make the limits diverge in the last 
four expressions. Thus, the m-functions have to be chosen such that a form similar to 
f ll04p . fll05p . fll06p and (I107P is recovered, which is achieved by choosing 

m + (E) = and nC(E) = ~. (108) 

b 2 CLi 

These relations are very important because they relate the m-functions to the explicit solu- 
tion allowing us to evaluate the spectral density function. 
Substituting ([96]) to (j99l) in fflOSl) . we get finally 

. _ GzUr[a,y(R)]-iGtUt[a,y(R)] hm 

[ } iG+UnaMR)}-G-Ur[a,y(R)Y 1 ' 
m - (E) = H+U 2 -[a,y(-R)}-tH-U 2 + [a,y(-R)] 

1 1 tH-U 2 + [a iy (-R)}-H+U 2 -[a } y(-R)Y { ) 



C. Spectral density and physical interpretation 

Once the Weyl-Titchmarsh functions are known, it is straightforward to compute the 
spectral density using f fl9|) . f|109p and (IllOp . It is a complicated functions of the interatomic 
distance R, the potential strength g and the electric field strength F. Because the potential 
is not bounded, it is expected that the spectrum will be a continuum on E G R. Bound 
states on the other hand should become resonances for which the width is related to the 
decay time of these quasibound states. To analyze these assumptions, the spectral density 
is plotted in the following and compared to the case without electric field. The position 
of poles of p in the complex energy plane is also investigated. Finally, the "plunging" of 
the ground state resonance in the Dirac sea is discussed. Notice here that this section is 
included to show some qualitative features of the model; a more rigorous treatment would 
be required to understand all of these results and is presently under study. 
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(a) 



(b) 







(c) 



(d) 



FIG. 5: Spectral density for a two Dirac delta potential for varying values of the electric 
field strength F. The value of the potential strength is fixed to g — 3.0, the value of the 
interatomic distance is R = 1.0 while the the electric field is (a) F = 0.2, (b) F = 0.5, (c) 

F = 1.0 and (d) F = 2.0. 



1. Dependence on field strength and interatomic distance 

In Figure [5j the spectral density for different values of the field strength is shown and 
compared to the case without electric field. As discussed previously, the system has two 
bound states when F = (if the delta potential strength g lies in the interval —4mc 3 R + 
2c\J Am 2 c 4: R 2 + 1 < g < Amc'R + 2cy/4m 2 c 4: R 2 + 1). As the electric field is turned on, the 
bound states become Stark-shifted unstable states: the ground state and first excited states 
are shifted in opposite direction while their width increases. This behavior is very similar 
to the non-relativistic results [32] . 
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It can be understood mathematically as the occurrence of resonances, that is poles on the 
unphysical sheet of the analytically continued resolvent operator 22J (or Green's function) 
given by (H — E)~ x (assuming of course that it admits an analytic continuation to the 
lower half complex plane in E). Recall that the resolvent operator has the same analytical 
structure in the complex energy E plane as the "trace spectral density" because the former 
is given by 49 ] 



ip + (x)ip (y) 



, for x < y 



g«(»,ioH "•;igr";if ) ' fOT ;:;. m 

m+(E)-m-(E)' iUI X ^ » 

where ^(x) = M ± (a;) + m (E)v (x). The solutions and v ± have no poles in E, so the 
analytical structure in the energy plane would come from the terms: 

1 and f(E) m -(_E) 



m+(E) -mr(E) m+(E)-mr(E) 

which are included in p. When there is no electric field, there are poles lying on E e M 
(bound states) and are solutions to fl56l) . When the electric field is nonzero, these poles 
move to the unphysical sheet and becomes eigensolutions to a non-hermitian operator 55] 
which allows complex valued eigenvalues. The effect of these poles shows up in the presence 
of peaks in the spectral density (if the poles are close enough to the real axis). 

The continua also show the main features of the non-relativistic results: as the field 
strength decreases, the oscillation increases and average out to the zero-field case in the 
limit F — y 0. For each bump of the continua corresponds a pole in the complex energy 
plane. It would be interesting to study if the branch cuts on (— oo, — mc 2 ] and [mc 2 , oo) 
when F = degenerates into these poles as the electric field is turned on. Physically, these 
resonances were interpreted in |32| as an example of Ramsauer-Townsend resonances. The 
figure also shows that there are two frequencies of oscillation. The first one corresponds 
to the backscattering of the electron by the two delta potential wells (inducing a smaller 
frequency of oscillation and also existing in the case without electric field) while the second 
one is associated to the scattering on the whole system, including the electric potential 
(inducing a larger frequency for small enough F) 32]. 

In Figure [6j the spectral density for various interatomic distances and fixed value of 
electric field strength (F = 0.2) is presented. It is well-known that as R — > oo in the 
zero-field case, the two bound states becomes quasi-degenerate [32|. This is because the 
two potential wells become independent as the typical tunneling probability between them 
-Ptunnei — > as R — > oo. In the presence of the electric field, the two bound states are 
shifted and do not become degenerate, as shown in the picture. As the interatomic distance 
is increased, their energy shift becomes more important (the symmetric ground state is 
pushed to lower energy while the anti-symmetric excited state is pulled to higher energy) 
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(a) 



(b) 





(c) 



(d) 



FIG. 6: Spectral density for a two Dirac delta potential and an electric field of strength 

F = 0.2, for varying values of the interatomic distance R. The value of the potential 
strength is fixed to g — 3.0 while the the interatomic distance is (a) R = 0.1, (b) R = 0.5, 

(c) R = 1.0 and (d) R = 2.0. 



and their widths decrease. The latter is due to a larger tunneling barrier for the electron as R 
becomes larg er. All of this is very similar to what one would expect from the non-relativistic 
analysis 32]. It also suggests a new process to probe the negative energy continuum: by 
taking R large enough for a given value of F, the real part of the ground state resonance 
could go into a region where it overlaps with negative energy states. We will see in Section 
IV C 21 however that this interpretation is slightly too simplistic as avoided crossings may 
occur depending on the value of g. 
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2. "Plunging" deeper in the Dirac sea 

For a given value of the electric strength, there is a critical value for the interatomic 
distance R for which the resonance energy -E^ ound = ^(E^ round ) + i^(E^ mund ) associated 
to the ground state (that is lim^o ^w>und = -Aground where -Ground is a solution to f )56|) ) 
approaches the negative energy continuum of the free operator, that is ^i(E^ round ) > —mc 2 . 
As the interatomic distance is increased further, one would expect that the resonance would 
go down below —mc 2 where most of the negative energy continuum states belongs. This 
phenomenon is hinted by the results obtained for the spectral density at large interatomic 
distance. A few examples of this are shown in Figure [7] where the spectral density for two 
electric strengths (F = 0.2 and F = 0.5) and two interatomic distances are presented. As 
seen in these pictures, the ground state resonance seems to cross the Ramsauer-Townsend 
resonances and seems to keep going to more negative energy and larger width as the inter- 
atomic distance is increased. 

However, a careful analysis of the position of poles in the complex energy plane reveals 
that it is otherwise. For a given value of g and as R increases, there are two possibilities: 
either the ground state resonance goes below —mc 2 monotonically or either it induces a 
series of avoided crossings (the crossing of resonances was discussed in without electric 
field and in [25] for the non-relativistic case). An example of this is shown in Figure [S] and 
[9] for two values of potential wells strength (g = 2.0 and g = 3.0) and where the the electric 
field strength is fixed to F = 0.2 while the interatomic distance is varied by increment of 
0.01. The position of poles in the complex plane is obtained by solving numerically the 
equation 

[iG+U+[a,y(R)\ - G~l7f[fl,y(i2)]] [iH;U+[a,y(-R)\ - H+U^{a,y(-R)]] 
- [iH;U+[a,y(-R)]- H+U2M-R)]] [iGtU?[a,y(R)}-G-U^[a,y(R)}] 

= 0, (113) 

which corresponds to the denominator of the spectral density. The figures show how the 
position of the poles varies as a function of R in the complex energy plane. As seen in 
these pictures, for g = 3.0, the first avoided crossing occurs when the real part of the bound 
state resonance approaches —mc 2 . At that point, its imaginary part changes rapidly to 
take the value that the nearby resonance had before the avoided crossing. In some sense, 
the two resonances exchange their positions in the complex plane, which is why they were 
indistinguishable in the peak structure behavior of the spectral density. For g = 2.0, there 
is no avoided crossings: the ground state resonance keeps going to lower energy as the 
interatomic distance is increased. 
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(c) (d) 

FIG. 7: Spectral density for a two Dirac delta potential. The value of the potential 
strength is fixed to g = 3.0 while the other parameters are (a) F = 0.2 and R = 4.0, (b) 
F = 0.2 and R = 5.0, (c) F = 0.5 and R = 2.5, (d) F = 0.5 and R = 3.0. There is a peak 
that goes into the negative energy continuum as R is increased, it is the fourth one from 
the left in (a) and (b) and the second one in (c) and (d). 



The physical implications of this phenomena (such as whether the system becomes over- 
critical |7| when the ground state resonance reaches $t(Eg VOUnd ) < — mc 2 ), its extension to 
a real 3-D system and a systematic study of the poles behavior for all parameters will be 
presented elsewhere. 
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FIG. 8: Real and imaginary part of the poles in the complex energy plane for different 
values of R and for F = 0.2. In (a) and (c), we have g = 2.0 while for (b) and (d), we have 
g = 3.0. The avoided crossings are surrounded by rectangles while the crossings are 

encircled. 



VI. CONCLUSION 



In this work, the spectral density for the 1-D Dirac equation was evaluated for three 
different cases using the WTK method. The free case was treated as a consistent test for 
the other systems studied. The main new feature in comparison to the non-relativistic 
case was the appearance of the negative energy continuum, as expected. The double delta 
potential wells potential was then treated. By varying the parameters of the model, it 
was possible to investigate the behavior of the symmetric and anti-symmetric states which 
were characterized by their eigenvalues positioned in the mass gap (although for sufficiently 
strong or weak potential wells, one of the state disappeared in the continuum). Finally, 
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FIG. 9: Position of the poles in the complex energy plane for different values of R. In (a), 
we have g = 2.0 while for (b), we have g = 3.0. The arrows point towards the direction of 
increasing interatomic distance. The rectangle surrounds the region where the avoided 

crossings occur. 



the last case considered was the same system immersed in an electric field. This allowed to 
investigate the appearance of Stark resonances and to characterize their behavior. It was 
shown that they behave in a similar fashion to the non-relativistic case, except when R or 
F is large enough in which case new relativistic effects were observed. 

The main new physical feature of our analysis is the fact that a resonance having the 
properties of the ground state resonance can move into the negative energy sea (of the 
free operator) for a sufficiently large interatomic distance and electric field strength. As the 
interatomic distance is increased, this happens either by simply moving in this energy region 
or through a sequence of avoided crossings. This phenomenon is clearly absent from the non- 
relativistic treatment. This suggests that the system becomes overcritical [3] in that case and 
thus, can produce particle-antiparticle pairs. Strictly speaking, as soon as the electric field 
is turned on, electron-positron pairs start to be produced from the Schwinger mechanism. 
When there are no potential wells, the probability of producing a pair is very small and 
so is the number of pairs produced, unless the field reaches a certain large value. Our 
analysis seems to imply that the presence of the diatomic molecule would allow decreasing 
the threshold for pair production by using the Stark effect. However, our work considered 
only the single particle Dirac equation: a definite conclusion can only be reached by a 
Quantum Field Theory treatment, which will be the topic of a future publication. Finally, 
one should also note that only a few qualitative features of the system were presented in this 
work. An exhaustive study of the resonance positions in the whole parameter space will be 
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the topic of future investigations. 

From the mathematical point of view, our work combined two main mathematical the- 
ories: the WTK method and the Colombeau's generalized function theory. The former is 
well-known and has been used to analyze many physically relevant equations and to obtain 
analytical solutions, as done in our analysis. To the best of our knowledge, it is the first time 
the latter is used to treat point interactions in the Dirac equation. Although the final result 
obtained is not completely new, we hope our analysis puts it on a more rigorous ground. 
For instance, the usual way of dealing with the electrostatic point interaction is to integrate 
the Dirac equation around the point interaction, that is on [— e, e] (for V(x) = —g8(x)) and 
then, make the following assumption 40, 5^]: 

,0+ x 

dxi/;(x)6(x) = - [^(0+) + V((T)] . (114) 

Recall that ip has a jump discontinuity at x = and thus, the latter implies certain assump- 
tions on the distribution representative of the potential and wave function. In our analysis, 
the Dirac equation is interpreted in the sense of Colombeau's distributions. We then use a 
symmetry argument and the general properties of self-adjoint extension of point interaction 
to fix the ambiguity in the product of distributions, allowing us to "derive" a relation similar 
to (I114p and to obtain boundary conditions for a self-adjoint operator consistent with an 
electrostatic potential. 



Appendix A: Colombeau's generalized theory for the treatment of point interactions 
in the 1-D Dirac equation 

The Dirac equation is given by 

Ei/j(x) = [ic<jy d x -a z mc 2 + V{x)]i){x). (Al) 
For the sake of this argument, we consider only one delta distribution interaction given by 



V(x) = —g5(x). 



(A2) 



The generalization to two or more delta interactions is straightforward. The goal is to 
find the boundary conditions for the wave function at x = 0. We do so by working in 
Colombeau's algebra which non-canonically embeds the set of Schwartz' distributions, al- 
lowing the multiplication of distributions, 45M48| . First, it was argued in 39j that contrary 
to the Schrodinger equation, where the wave function is continuous (its derivative has a 
jump discontinuity), the wave function computed with the Dirac equation has a jump dis- 
continuity. Therefore, when considering the term V(x)ip(x), we are dealing with a product 
of distributions. 
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The first step to find the boundary conditions is to rewrite ( lAlj) componentwise. This 
yields 



Eipi(x) = cd x ip 2 (x) — mc 2 ipi(x) + V \x)ipi{x) , (A3) 
Eip 2 (x) = -cd x ipi(x) + mc 2 ip 2 (x) + V{x)^ 2 (x), (A4) 

which are true in the sense of distributions. Then, we are seeking discontinuous solutions of 
the form 47] 



ipi{x) = Aip 1 (x)H(x) + ip 1 (x) where Aipi(x) = ipi(x) — ip^ (x), (A5) 
ip 2 (x) — Aip 2 (x)K(x) + ip 2 (x) where Atp 2 (x) = iptix) — ^{^i (A6) 

where ip 12 are smooth continuous functions defined in R and solution of the free (g = 0) 
Dirac equation, H and K are Heaviside generalized functions. A priori, the latter are equal 
in the weak sense only (K « H): that is there exists (i^ e ) < e <i, (H e ) 0<e< i, two distinct 
representatives respectively of K and H, then for all G P(R) (set of C°°(R)-functions 
with compact support): 

lim / K e (x)cj)(x)dx = K(<j>) = / 4>{x) = lim / H e (x)<f>(x)dx = H(4>) 

e ^°jR Jr+ e ^°jR 

such that V0 6 2?(R) : 

lim / (fC(x) - F e (a;)U(a;)da; = 0. 



In the strong sense, these last equalities would require K t {x) = H £ (x) (see [47j, for details 



on these definitions). In the strong sense 47j, we assume that H ^ K. We will show a 



posteriori that they are also equal in the strong sense. Substituting (1A5|) and (1A6|) in (1A3|) 
and (IA40 . we get, in Colombeau's distributional sense 

A[$,E\ := -EAfa(x)H(x) + c[d x Aifa(x)]K(x) 
+cAip 2 (x)K'(x) - mc 2 A'ip 1 (x)H(x) 

-gL'(x)A^j 1 (x)H(x) - gL'(x)^(x) = 0, (A7) 
B[$,E\ := -EA^ 2 (x)K(x) - c[d x Afa(x)]H(x) 
-cA^ 1 (x)H'(x) + mc 2 Aijj 2 {x)K{x) 

-gL'(x)Aij; 2 (x)K{x) - gL'(x)ip 2 (x) = 0, (A8) 

where L is also a Heaviside function and H' = d x H, K' = d x K, L' = d x L are delta functions. 
The fact that ip 12 are solutions of the free Dirac equation was also used to cancel some 
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terms. Then, to find the boundary conditions, we apply G Z?(R) to (1A7j) and (1A8|) : 

v4 e [^,-E] := / Aip 1 (x)H e (x)(f)(x)dx + c / [^A^^)]^^)^^)^ 
Jr Jr 

+c / A-0 2 (x)K' e (x)(j)(x)dx — mc 2 / Ai/j 1 (x)H e (x)(f)(x)dx 
Jr Jr 

—g / L' e (x)Aipi(x)H £ (x)(f)(x)dx — g / L' e (x)^jf(x)(/)(a;)(ia; 
= 0, (A9) 

Be[ip,E) := -E / Aip 2 (x)K e (x)(f)(x)dx - c / [9 I A^i(x)]F e (x)0(a;)dx 
Jr Jr 

—c / Avpi(x)H' e (x)(J)(x)dx + mc 2 / A^ 2 (a;)-ft' e (x)0(x)(ix 
Jr Jr 

—g / L' e (x)Aip 2 (x)K e (x)(f)(x)dx — g / L' £ (x)ip2(x)(J)(x)dx 
Jr Jr 

= 0. (A10) 

Performing the integration and taking the limit e — >■ 0, we get 

0(0) [c^+(0) - c^j(0) - Wi + (0) - <?(1 - a)Vr(O)] 

[(£ + mc 2 ) A^i(x) - cc^A^Or)] 0(x)dx = 0, (All) 



0(0) [-c^ + (0) + c^r(0) - ^+(0) - (7(1 - 6)V>2 (0)] 

[(£ - mc 2 )A^ 2 (x) + c^A^i(x)] 0(x)dx = 0, (A12) 



where we used the fact that ■0 ± are smooth and where a, b are real constants. These constants 
are obtained from Colombeau's theory which states that for all G Z>(R) 

lim / L' e (x)(p(x)H e (x)dx = a L' e (x)(p(x)dx = a0(O), (A13) 
^°Jr Jr 

lim / L' e {x)<j)(x)K e (x)dx = b L' e {x)<f){x)dx = 60(0), (A14) 
Jr Jr 

for real constants a, b which value depend on the local structure of distributions. Now, 
we can get rid of the integral terms in (lAllj) . (IA12I) by using again the fact that tpf 2 are 
solutions of the free Dirac equation and that the equality is satisfied for all 0. We get: 

c^ + (0) - c^(0) - ga^+(0) - g{l - a)^(0) = 0, (A15) 
-q#(0) + c^r(0) - ^+(0) - g(l - 6)VJ(0) = 0. (A16) 

Thus, we have the boundary conditions or jumps conditions on the wave function since we 
have lim^o ip(±e) = , ± (O) from (IA5I) and (IA6I) . They are ambiguous because the value of 
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a, b is unknown: they should be fixed by some other physical or mathematical arguments. 
This is done in the following. 

There are many ways to implement point interactions in the 1-D Dirac equation and each 
corresponds to a different self-adjoint extension of the Dirac operator. It was shown in 41] 
that the most general self-adjoint extension is given by a four-parameter family (a, 5, 7, (3): 
for each quadruple corresponds a different physical system. This is believed to be the 
main reason for the ambiguities arising in the definition of the boundary conditions. More 
precisely, the operator implementing point interactions is given by the free Dirac operator 
H on the domain 



4l| 



V{H) = G {^(RUO}) n AC loc (l\{0})} ® CV(0 + ) = A^((T)} (A17) 

where W 2,1 is the Sobolev space, AC\ QC is the set of absolute continuous function and A is a 
2x2 matrix implementing the following transformation: 



<Mo+) =u [<tyi(o-) + 7^2(0-)] . 

^ 2 (0 + )=w[/3Vi(0-)+a^(0-)] 



(A18) 
(A19) 



where a, 8, 7, j5 £ R, a8 — = 1, co G C and \u\ = 1. This choice of parameters guaran- 
tees that the Dirac operator is self-adjoint and reduces to the usual delta interaction of the 
Schrodinger equation in the non-relativistic limit. Therefore, the jump condition derived 
previously in (IA15j) and flA15j) has to be consistent with these conditions to keep the self- 
adjointness of the operator. Moreover, it should also fulfill an important "physical" condi- 
tion: the charge conjugation invariance. The latter guarantees that the particle-antiparticle 
duality is well-defined. Both of these requirements allow to fix the value of a, b uniquely, as 
shown in the following. 

First, we show that K = H in the strong sense by demonstrating that a = b using 
the charge conjugation invariance. It can be easily found from the time-dependent Dirac 
equation in our representation that the wave function transforms as Ctp(x) = a x ip*(x) under 
the charge conjugation operator C. Recall that the potential then transforms as g — > —g. 
Then, (I A 1 5 [) and flA16j) are written in the same form as (IA18j) and flA19j) . which yields 



Vi(0" 



1 + 


abg 2 


-1 




c 2 _ 




( 


1 + 


abg 2 ' 


-1 




c 2 _ 




( 



\- b -^)Mo-)- 9 -Mo-) 
i-f*?^W) + f*(o-) 



(A20) 
(A21) 
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and 



5 = 


q 2 ab 
c 2 


-l - 


a = 


q 2 ab 

l + — 

c 2 


-1 r 



1 - 



g%{l-a)- 



g 2 a{\ - b) 



7 

CO 



1 + 



g 2 ab 



-1 



9 
c 



-0, 



1. 



(A22) 

(A23) 

(A24) 
(A25) 



We require that the operator defined on the domain T>(H) is invariant under the charge 
conjugation symmetry, which implies that the matrices A c = CAC -1 should be the same as 



A 



such that the following equality holds: 

5 7 



A := u 



(3 a 



A 



5 —7 
—j3 a 



(J, 



(A26) 



Note that the sign of /3,7 in A c changes because g — > —g. The last equation implies that 

(A27) 





5 7 


* 


a —(3 




(3 a 


= u 


—7 5 



which is fulfilled if and only if a = b. 
From this last argument, we obtain 



1 + 



2 2 1 ~ 1 
9 a 21 



1 - 



g 2 a{\ 



a, 



7 
to 



1 + 



g 2 a 2 



-1 



g 

c 



-p, 



1. 



(A28) 

(A29) 
(A30) 



which is still ambiguous because the value of a is unknown. However, recall that the self- 
adjoint extension is defined such that the condition a5 — 7/? = 1 is fulfilled. Substituting the 
last three equations in this condition, we find that a is the solution of the cubic equation 



1 2 C 2 C 2 



> > 0, (A31) 
2 g 2 2g 2 

which has two complex solutions and one real solution given by According to 

Colombeau's theory, the constant value should be a real number so we select a 
fixes the value of a, b uniquely and unambiguously. 

To summarize, we have the following boundary conditions: 

,2 



\. This 



Vi(0" 
^2(0- 







-1 




1 + 


g 2 ' 




( 


4c 2 _ 








/ 


-1 




1 + 






( 


4c 2 _ 







1 - 



9 
Ac 2 
2 



^i(o-) - ^2(0-) 

c 



i-^)^(o-) + f^i(o- 



(A32) 
(A33) 



33 



and also, it is convenient to invert these to get 

^2(0-) = 



9 2 ' 


-1 




4c 2 _ 




V 4c 2 




-1 




9 2 ' 






4c 2 _ 




A 4c 2 



^i(0 + ) + ^ 2 (0 + ) 
c 

c 



(A34) 
(A35) 



These boundary conditions are consistent with the ones found in 57|, |59) for static electric 
potential which were derived by using other but similar techniques. 

We conclude this Appendix by a remark about the definition of the product of distri- 
butions. In the framework of nonlinear hyperbolic systems, Dal Maso, LeFloch and Murat 
have proposed a definition of such products (typically H x 5) as a bounded Borel measure 
dependent of a path $ (which has several properties 60[]) as follows. Assuming that if) is 
regular expected at xo where it has a jump, and assuming F to be a regular function, then 
the product F(ip)ifj' is defined as the following measure. If ip is continuous on a Borel set 
B, then 



F{$)%l/(x)dx 



B 



At Xo, where ip is discontinuous 

[F(^'Ux ) = J F($( S ;<Kx ),^(a; +))— (s;^)^(x+))ds 
where by definition $(0; iP(xq ), iP{xq )) = 4>(xq) and $(l; %Jj(xq), tp(xQ )) = ^(^o)- ^ 

F($(s;^o, Vi))-^-(s; ipo,ipi)ds x <J Xo 



Now in the case, where F is the identity, we get 

w\* -- 



1 5$ 



9s 



1 



That is: "if 5 = S/2" , which is the same result obtained from Colombeau's theory. Details 
can be found in 



60|], [6lJ. 
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